Transient dynamics of a single molecular transistor in the presence of local electron–phonon and electron–electron interactions and quantum dissipation

We consider a single molecular transistor in which a quantum dot with local electron–electron and electron–phonon interactions is coupled to two metallic leads, one of which acts like a source and the other like a drain. The system is modeled by the Anderson-Holstein (AH) model. The quantum dot is mounted on a substrate that acts as a heat bath. Its phonons interact with the quantum dot phonons by the Caldeira–Leggett interaction giving rise to dissipation in the dynamics of the quantum dot system. A simple canonical transformation exactly treats the interaction of the quantum dot phonons with the substrate phonons. The electron–phonon interaction of the quantum dot is eliminated by the celebrated Lang-Firsov transformation. The time-dependent current is finally calculated by the Keldysh Green function technique with various types of bias. The transient-time phase diagram is analysed as a function of the system parameters to explore regions that can be used for fast switching in devices like nanomolecular switches.

With the advances in fabrication techniques in nanotechnology, together with the availability of ultrafast experiments to analyze the out-of-equilibrium quantum dynamics of many body interacting systems, research in molecular electronics has gathered unprecedented momentum in recent years. Indeed, ultrafast experiments have been employed to study the charge transport in a three-terminal device like a single molecular transistor (SMT) [1][2][3][4][5][6][7] . Recently, Cocker et al. 8 have studied the charge transport dynamics through a singly-occupied localized orbital, like highest occupied molecular orbital (HOMO) in pentacene molecule using the terahertz-scanning tunneling microscope (THz-STM) technique. Many experimental groups have fabricated molecular devices using the C 60 molecule [9][10][11] and found that the current is controlled effectively by tuning the gate voltage. Also, many researchers have reported the non-equilibrium effects of electron-phonon (e-p) interaction during the charge tunneling in molecular devices like phonon-assisted tunneling transport [12][13][14][15][16][17][18][19][20][21][22] , hysteresis-induced bistability 23 , local heating 24 , molecular switching 25,26 , and negative differential conductance [27][28][29] . Chen et al. 30 have studied the polaronic effects on the current density using the Keldysh formalism. They have reported that the current density decreases with increasing the e-p interaction and also the formation of phonon side peaks in the spectral density. Later, Raju and Chatterjee 31 have studied the effects of quantum dissipation due to substrate in an SMT device with electron-electron (e-e) and e-p interactions and reported that the dissipation enhances the current density. In recent work, we have investigated the magneto-transport properties 32 of an SMT system and explained that it can be used as a spin-filtering device at zero temperature. We have also shown the effect of e-p interaction on the transport properties at finite temperature in 33 .
Of late, the transient dynamics of mesoscopic systems has attracted considerable attention. In this context, dissipative driven quantum rings, dynamical Franz-Keldysh effect, molecular electronics have been extensively studied. The theoretical investigation of transient dynamics in molecular electronic devices with strong e-p interaction is a challenging task. It is essential from the point of view of the application of these devices for quantum computing. Jauho et al. 34 have analyzed the transient dynamics through a non-interacting quantum dot (QD) response to harmonic and sharp square-shaped voltage pulses by using the Keldysh Green function method. Schmidt et al. 35 have investigated the time-dependent transport through QD with weak Coulomb interaction using the Anderson model within the framework of a mean-field approximation and used the Monte Carlo technique to calculate the current density. Many researchers [36][37][38] have investigated the transient dynamics of QD modeled by the impurity Anderson model in the Kondo regime by employing different approaches like the time-dependent noncrossing approximation method 39,40 , time-dependent density matrix renormalization group technique 41 , functional renormalization group approach 42 . In reference 43 , they studied the transient dynamics of a single molecular junctions modeled with the Anderson-Holstein model using the auxiliary-mode expansion method . In the present work, we study the transient dynamics of an SMT device with e-e and e-p interactions and quantum dissipation by employing Keldysh formalism. Here we explore the transient dynamics of the SMT device with two different time modulations, namely, the harmonic pulse modulation and the upward pulse modulation. Figure 1 represents a block diagram of an SMT device with a molecule or QD is coupled to two metallic leads, a source ( S ) on the left and a drain ( D ) on the right. The entire structure is mounted on an insulator (substrate), which acts like a phonon bath. The electron energy levels in the QD system can be tunable by changing the gate voltage V g . The system is driven by the bias voltage V b which is considered to be time-dependent. To model the system under consideration, we consider the extended AH model Hamiltonian by incorporating the Caldeira-Leggett 44 term to take care of the dissipation. The Hamiltonian H can be written as:

Model
Here H l is the Hamiltonian describing the leads and is given by where n αkσ (= c † αkσ c αkσ ) gives the number of the conduction electrons in the lead α with momentum k , with time-dependent energy, ε αkσ (t) = ε 0 αkσ + � α (t) , � α (t) being the time-dependent bias which is equal to V b (t) and c αkσ (c † αkσ ) is the corresponding electron annihilation (creation) operator. The second term in the Hamiltonian H c models the central region or the QD and is given by where n dσ = c † dσ c dσ describes the number operator for the central dot electrons with time-dependent energy ε d ,c dσ (c † dσ ) represents the dot electron annihilation (creation) operator, U is the onsite e-e interaction, b † (b) is the creation (destruction) operator for a QD phonon with frequency ω 0 and is the measure of e-p interaction strength. The third term in Hamiltonian, H t represents the tunneling coupling between the QD and leads with coupling strength V k (t) , can be written as Finally the H B term describes the substrate phonons and their interaction with the phonon of the QD, is given by (2) www.nature.com/scientificreports/ where x j and x 0 denotes the generalized coordinate of the substrate and dot oscillator, ω j gives the frequency of the j th substrate oscillator and β j represents the linear coupling strength between the j th oscillator of the substrate and dot oscillator.

Formulation Elimination of the interaction with bath phonons. The linear interaction between the dot and bath oscillators can be decoupled by performing the following unitary transformation
As a result, the local phonon frequency of QD ω 0 gets renormalize to ∼ ω 0 , is given by is the change in ω 2 0 due to the interaction with the bath phonons. The decoupled Hamiltonian is given by Here onwards we will discard the bath Hamiltonian because it merely adds constant energy to the system. We considered the spectral function J(ω) to describe the dynamics of the substrate phonons which is given by For all frequencies ω , the spectral function obeys the relation in the strict Ohmic case: J(ω) = 2m 0 γ ω , where the ohmic damping coefficient: The pure Ohmic spectral density calculated above, however, is not very practical because it diverges in the limit of ω → ∞ . In order to salvage the situation, we have considered the Lorentz-Drude form 45 for J(ω) with the cut of frequency ω c , in the limit ω → ∞ , J(ω) approaches the value zero and in the small-ω limit, we can obtain the pure Ohmic spectral density. The final expression for  www.nature.com/scientificreports/ where is the Fermi distribution of the lead α and µ α is the corresponding chemical potential which is related to V b as: µ α = eV b δ αS , here δ αS is the Kronecker delta function and Ŵ α is given by with Within the local polaron approximation we can replace, represents the phonon population, and at zero temperature Ṽ k = V k (t)e − 2 /2 . Here, our numerical computation is done at zero temperature which also evident that we have considered that the local phonons are always in thermal equilibrium with the bath phonons. Where ρ S(D) represents the constant density of states in the source (drain). Within the wide-band approximation, ρ s,D (ε) are independent of energy, and N d,σ (t) is the occupancy of the energy levels in QD, which is given by is the spectral function of the dot, which is calculated by using the Keldysh Green function method as where G r d,σ (t, t 1 ) is the retarded Green function given by Using the equation of motion method, we get the expression for G r (t, t 1 ) as To get the above expression we have treated the e-e interaction in the mean field level, i.e.
Here we investigate the time-dependent current through the system for two time modulations.

(i) Harmonic pulse time modulation
In the case of harmonic pulse modulation, � α (t) = � α cosωt and substituting in Eq. (18) and after some algebraic manipulation, we obtain the imaginary part of corresponding spectral function as, where J k (� α /ω) is the Bessel function of the first kind and ω = 2π/T,T = (πℏ/ Ŵ) being the time period of the harmonic bias voltage. (

ii) Upward pulse time modulation
In this case, we consider � α (t) = � α u(t) , where u(t) is the Heaviside function and the corresponding spectral function is calculated as Substituting the imaginary part of the above function in Eq. (13) finally, yields the time-dependent current density.

Results and discussions
In the present work, we consider a single energy level in QD with energy ε d = 0.5, and symmetric coupling (Ŵ S = Ŵ D = Ŵ/2 ) of the dot level with the source and drain. Here, we measure all the energy quantities in terms of dot phonon energy ℏω 0 , time in units of ℏ/ Ŵ and in rest of the paper we set ℏ = 1 . Here we have studied the normalized current J/J 0 = (J S − J D ) with J 0 = eŴ/2ℏ for two different input pulses. Here we have taken eV g = 0.3, γ = 0.03, ω c = 5, Ŵ = 0.2, eV b = 0.5, k B T = 0, U = 3 , D = 0 and S = 5 . We first present our results for the Harmonic pulse and then we will discuss the results for the upward pulse.
Harmonic Pulse modulation. Figure 2 represents the behavior of occupancy of the dot as a function of gate voltage for various e-p coupling strengths at t = 10 . We have calculated the values of N d↑ and N d↓ self-consistently subjected to the condition that N d↑ + N d↓ = 1. The occupancy clearly displays the staircase structure due to phonon side bands. The up-spin occupancy N d↑ increases with increasing gate voltage while N d↓ shows a decreasing behavior. In Fig. 3, we have plotted the normalized current density versus time for different values with J 0 = eŴ/2 . Here, one can observe that the oscillations in the normalized current decrease with increasing time and reach a steady state at a particular time called transient time. We can observe that the transient time decreases with increasing e-p interaction strength and it can be admitted to variations in the hopping parameters. The reason can be explained as when an electron travels from the source to the dot, and it interacts with the phonons as a result,     www.nature.com/scientificreports/ it forms a polaron. The polaronic effects can be clearly observed from Eqs. (12) and (15). The current density approaches the steady-state faster as the e-p interaction strength increases as a result the transient time decreases. In Fig. 4, we have plotted the current density as a function of time for different values of the damping rate. One can observe from the inset that the damping rate slightly increases the current density, which is expected behavior. We have presented the current density versus time result for different gate voltages in Fig. 5 to see the    Fig. 6. It is evident that as time starts, the current density fluctuates much more with increasing gate voltage, while as time advances sufficiently, the fluctuations subside. In Fig. 7, we show the current density versus result for different time values. One would expect that the normalized current density would decrease with increasing due to polaronic effects. However, many fluctuations seem to occur at a small time and one has to actually wait for some time to get the stable value of the current density. In Fig. 8, we show the three-dimensional plot and the contour map of the current density as a function of t and . It is clear from the figure that reduces the transient time, same as observed in Fig. 3. Figure 9 represents the contour map of the current density as a function of e-e interaction and time. It is evident that the current density decreases with increasing e-e interaction at the smaller time scale, but after the current reaches its steady state, the effect of e-e interaction is very small.

Upward Pulse modulation.
We now present the results for the upward pulse modulation for which we set the values of the system parameters as:Ŵ = 0.01, eV g = 0.6, eV b = 0.5, k B T = 0, ℏω 0 = 1 , D = 0 and S = 5 . Figure 10 shows the effects of e-p coupling on the dot occupancy as a function of gate voltage. Again the staircase structure is visible, though the behavior of N d↑ and N d↓ are opposite as a function of tunable parameter eV g . However, the e-p coupling strength reduces the occupancy for both up and down spin electrons to maintain the total occupancy is equal to one. We show the current density behavior as a function of time with different e-p coupling values in Fig. 11. The e-p interaction is seen to reduce the amplitude current oscillations and transient time as observed in the case of harmonic pulse modulation.
In Fig. 12, we plot the current density as a function of time with different gate voltages and at a particular value of e-p coupling. Comparison of the transient time for several values of the gate voltage shows that it is smaller at higher values of the gate voltage. Thus the gate voltage also reduces the transient time. For a better understanding in Fig. 13, we show the current density versus gate voltage result for different t values. Clearly, the J/J 0 reduces    www.nature.com/scientificreports/ with staircase structure with increasing gate voltage because as the gate voltage increases, the dot level moves out of the bias range. Also, the width of the plateau region with constant current is increasing for higher values of time. In the plateau region, the device can be used for the switching applications. We plot the current density as a function of at different times in Fig. 14. It is interesting to see that the oscillations in the J/J 0 reducing with increasing the e-p interaction strength and at the same time, current density decreasing with . Figure 15 represents the three-dimensional and contour plot of current density as a function of time and . Thus, one can conclude that reduces the transient time and also as increases, the current density saturates when t is beyond 50. Thus, one can conclude that reduces the transient time and current becomes almost constant at higher values of due to a strong polaronic effect. In Fig. 16, we have shown the contour map of the current density as a function of e-e interaction and the time to see the parameter region to obtain the steady state current.

Conclusions
In this work, we have studied the transient dynamics of an SMT device with e-e and e-p interactions and quantum dissipation for the harmonic pulse and upward pulse modulations. The system is modeled by the Anderson-Holstein-Caldeira-Legget Hamiltonian. The linear interaction of dot phonons with the substrate phonons, which introduces the dissipation, is treated with an exact transformation that essentially modifies the QD phonon frequency. Next, the e-p interaction is decoupled from the Hamiltonian by applying the Lang-Firsov unitary transformation and finally, we employed Keldysh formalism to calculate the current density expression. We investigated the variation of the transient time with respect to the system parameters like e-p interaction strength, damping coefficient and eV g for both the harmonic and upward pulse modulations. The most important result that we have perceived is that the transient time decreases with increasing the e-p and e-e interactions. Also, we have shown the parameter region with constant current, where the device can be used for switching applications. The current density decreases with increasing e-p interaction strength due to the polaronic effect. Our theoretical results can be easily tested with the experiments on a single molecular transistor, nanomolecular switching devices.